Bridging the Time Scale Gap with Transition Path Sampling

Transition path sampling is a methodology which overcomes both the long timescale problem and the lack of prior knowledge about transition mechanisms. Here we briefly review the basic principles of transition path sampling, illustrate its application using autoionization in liquid water, and emphasize the capabilities and limitations of the methodology. 1 Why transition path sampling is needed Many interesting processes in nature are characterized by the presence of different relevant time scales. In a chemical reaction, for instance, the reaction time can be many orders of magnitudes longer than the molecular vibration period usually measured in units of femtoseconds [1]. Such a separation of time scales creates serious problems for the computer simulation: on one hand the resolution in time needs to be fine enough to capture the properties of fast motions (such as molecular oscillations) and on the other hand the simulation must be extended to times longer then the longest relevant time scale in order to observe the events of interest (such as chemical reactions). This is the notorious time scale gap problem addressed in this conference. It is not only a problem in chemical physics. For example some comets exhibit rapid transitions between heliocentric orbits inside and outside the orbit of Jupiter [2]. While the transition, during which the comet transiently orbits Jupiter for a few periods, is swift, many revolutions of the comet around the sun can occur between transition. Often, widely separated time scales are caused by energy (or free energy) barriers preventing the system from quickly visiting a representative sample of pertinent configurations. In the past, many efficient computer simulation techniques [3] such as umbrella sampling [4], the multiple histogram method [5], and, most recently, the Laio-Parrinello approach [6], to mention just a few, have been developed to overcome such barriers and sample the free energy surface for a specified control parameter (or order parameter). While biasing schemes of this sort can be used to determine structural quantities such as chemical potentials and equilibrium constants, they are of limited use if one wishes to study the 2 Christoph Dellago and David Chandler dynamics or kinetics of rare transitions. In this case, detailed knowledge of the underlying mechanism in terms of a reaction coordinate is necessary to apply standard techniques such as the reactive flux method [7–9]. For numerous problems of current interest, however, reaction coordinates are not known, and in applying the techniques mentioned above one has to rely on guessing the relevant, possibly collective, variables. If the selected degrees of freedom do not capture the essence of the mechanisms under study, such biasing schemes are bound to fail. Fortunately, there are diagnostic tools, such as committor distributions (to be described later in the article), which can be used to detect if the relevant degrees of freedom have been identified correctly or not [10,11]. But what can one do when such failure has become apparent? To answer this question, it is important to distinguish whether the system under consideration is simple or complex. The former is where dynamics is dominated by energetic (as opposed to entropic) barriers, and the topography of the energy landscape is not excessively complicated. In this case, insight can be obtained by locating minima (stable states) and saddle points (transition states) on the potential energy surface. This basic idea has been successfully implemented in various methods including eigen vector following [12,13], the nudged elastic band method [14], and hyperdynamics [15]. In complex (i.e., non-simple) systems, however, the potential surface can exhibit a huge number of distinguishing features, such as local minima, maxima and saddle points. Explicitly enumerating all these features is impractical. Further, dynamical bottlenecks need not coincide with these features in any straightforward way. This fact is illustrated by the ”golf course” landscape, in which an entropic barrier hinders the system from finding its energetically most favorable configuration [16]. The failure of methods relying on searching the potential energy surface for stationary points can be illustrated even more drastically in entropy driven phase transitions, such as the freezing of hard spheres experimentally observed in colloidal suspensions [17,18]. In such processes all configurations accessible to the system are isoenergetic and any transition is purely driven by entropic imbalances. No potential energy minima or saddle points exist, but there still are stable states separated by a free energy barrier and rare transitions between these stable states can occur. So what then can we do for a complex system where we don’t know the mechanism and cannot find it by searching for specific points on the potential energy surface? Transition path sampling offers an answer by considering trajectories instead of single configurations [19–21]. This change in perspective permits the generation of rare transitions between stable states without prior knowledge of mechanisms or reaction coordinates. Rather than requiring such prior information as an input, transition path sampling can help in finding it. In the following section we briefly review the essential ideas of transition path sampling. More detailed descriptions of the method including several illustrative examples are given in two recently published review articles [10,11]. As with any newly developed technique, there are misconceptions about what transition path sampling can and cannot do. As such, we follow our summary of the basic principles of Transition Path Sampling 3 transition path sampling with a discussion of its capabilities and its limitations, with emphasis on how these limitations might be surmounted in the future. 2 How transition path sampling works 2.1 Probabilities of trajectories The transition path sampling method is based on a statistical mechanics of trajectories in which every trajectory x(t) of length t is assigned a statistical weight P [x(t)]. The set of all pathways x(t) consistent with the path probability P [x(t)] is called the transition path ensemble. Here x is a possibly high dimensional vector including all variables necessary to specify the state of the system under study. For instance, in a molecular system x may consist of coordinates and momenta of all particles. The trajectory x(t) is a sequence of such states, in which x0 denotes the first state on the trajectory and xt denotes the last one. For practical reasons x(t) is represented by a or chain of states, but in principle x(t) can be thought of as a general trajectory, which can be continuous or discrete depending on the process one intends to study. The form of the probability functional P [x(t)] specifying the weight of a given path x(t) in the transition path ensemble depends on dynamical rules governing the time evolution of the system. Let’s for simplicity assume that the system evolves according to some set of equations of motion (for instance Newton’s equations of motion) and that the dynamics has the Markov property, i.e. the state x of the system at time t completely determines the probability to find the system in a certain other state x′ a short time later. Then, the probability of a particular trajectory to be observed is: P [x(t)] = %(x0) L−1 ∏ i=0 p(xi∆t → x(i+1)∆t) , (1) where we have imagined that the pathway is represented by an ordered chain of L states and p(xi∆t → x(i+1)∆t) is the (properly normalized) conditional probability to observe the system in state x(i+1)∆t at time t + ∆t provided it was in state xi∆t a time ∆t earlier. In the case of Newton’s equations of motion (and of any other deterministic set of equations of motion) the path probability P [x(t)] consists of a product of Dirac delta functions describing deterministic trajectories flowing from the respective initial conditions. If the system under study is more conveniently described by a stochastic equation of motion (such as Langevin’s), the short time transition probability p(xi∆t → x(i+1)∆t) is not singular and more than one trajectory is allowed to emanate from the same initial condition x0. In any case, the path functional describes trajectories generated with a particular dynamical rule. Therefore, members of the transition path ensemble are physical trajectories (as opposed, for instance, to artificial minimum energy pathways) that can be used to study kinetics and dynamics. In (1), %(x0) denotes the distribution of starting points. This distribution of initial conditions might be the canonical distribution if one is studying a 4 Christoph Dellago and David Chandler system in thermal equilibrium. In other situations, %(x0) may represent a specific non-equilibrium distribution of initial conditions generated in a particular experiment. 2.2 Defining the transition path ensemble Let us now restrict the transition path ensemble to trajectories starting in a certain region A in phase (or configuration) space and ending in a different region B: P [x(t)] = hA(x0)%(x0) [ L−1 ∏ i=0 p(xi∆t → x(i+1)∆t) ]